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Abstract 


A recursive technique for modeling and estimating a two-dimensional signal 
contaminated by noise is presented. A two-dimensional signal is assumed to be an 
undistorted picture, where the noise introduces the distortion. Both the signal and 
the noise are assumed to be wide-sense stationary processes with known statistics. 
Thus, to estimate the two-dimensional signal is to enhance the picture. 

The picture representing the two-dimensional signal is converted to one dimen- 
sion by scanning the image horizontally one line at a time. The scanner output 
becomes a nonstationary random process due to the periodic nature of the scanner 
operation. Procedures to obtain a dynamical model corresponding to the auto- 
correlation function of the scanner output are derived. Utilizing the model, a 
discrete Kalman estimator is designed to enhance the image. 


vi 


JPL TECHNICAL REPORT 32-1596 



Two-Dimensional Signal Processing With 
Application to Image Restoration 


I. Introduction 

In theory, image enhancement utilizing classical filter- 
ing techniques does not seem to be difficult, since the 
processing of signals in one dimension can usually be 
extended to a two-dimensional case (using the concept 
of linear system theory). However, the application is 
cumbersome and often becomes impractical when a large 
amount of data is to be processed. Thus, it seems desir- 
able in many situations that a recursive technique of 
image enhancement be developed. 

In this report, we shall consider those images that are 
distorted by random noise and are best characterized by 
statistical procedures, such as specifying their first two 
statistical moments associated with the random process 
representing the brightness level. Thus, in this case, the 
image enhancement becomes a problem of statistical esti- 
mation and filtering, and to enhance a distorted image is 
to estimate the image. The input to the estimator is the 
output of a horizontal line scanner scanning the image 
with uniform speed. The procedure is first to develop a 
dynamical model whose response characteristic matches 


that of the scanner output in the statistical sense. The 
image is two-dimensional, while the scanner output is 
one-dimensional; thus, the model must exhibit the vertical 
correlation of the image. This correlation will be revealed 
by the oscillatory nature of the model responses. 

It is assumed that the image is characterized by a sta- 
tionary, two-dimensional autocorrelation function. How- 
ever, the scanner output, denoted by s(t), is nonstationary 
and nonseparable (Ref. 1) because of the scanners periodic 
movement. Consequently, no finite dimensional linear 
dynamic model representing s(t) exists. To remedy the 
nonseparability of the nonstationary process. s(t) 7 which 
is a very undesirable situation (Ref, 1), we shall generate 
a stochastic process whose autocorrelation function is 
stationary and which approximates the autocorrelation 
function of s(t ). 

Since we shall be dealing with the question of reali- 
zation of autocorrelation functions and thus spectral fac- 
torization, a brief background of spectral factorization 
is presented. 
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A. Spectral Factorization 

The concept of spectral factorization has become in- 
creasingly more important since Weiners original work 
(Ref. 2) on the subject Basically, spectral factorization 
determines the equations that describe a linear system 
when the system is driven by white noise and the covari- 
ance of the output is known. Whenever the covariance 
function of a process is driven by white noise via a sys- 
tem of differential equations of first order, we refer to this 
system as a dynamical model. More specifically, given a 
covariance function R(t, r), where t — and r^n for 
some fixed fi and n, the factorization problem is to deter- 
mine a realizable linear filter (differential equation model) 
that, when driven by white noise, yields R(f, r) as its out- 
put covariance. 

It is well known (Ref. 1) that, in general, no such reali- 
zation may exist. However, if its existence were guaran- 
teed, the representation (in some sense) would be unique 
(Ref. 3). In its most popular form, the spectral factoriza- 
tion would be confined to stationary situations. Then the 
corresponding dynamical model under consideration 
would be time-invariant, and the white noise forcing 
function must have started infinitely in the past (Ref. 3). 
This dynamical model would be asymptotically stable 
(Ref. 3). It is also desirable to deal with finite-dimensional 
dynamical models, implying that each linear model must 
possess a rational bilateral Laplace Transform. We can 
summarize the above discussion by the statement of 
Theorem 1, which we shall not prove (Ref. 4). 


with degree less than or equal to n — 1 and all roots in 
the left half of the s -plane, where and /Si are the real 

coeffi dents. That is, H(s) has all of its poles and zeros in 
the left half of the s-plane, 

B. Determination of the Output Covariance From a 
Linear Dynamical Model 

Consider the following dynamical model, given by 

x — A(t)x(t) + B(t)u(t) 
y(t) = C(t)x(t) (2) 

where x(t) is an n X 1 vector, u is an m X 1 vector, y is a 
scalar, A, B, and C are matrices of appropriate dimen- 
sions (not necessarily time-invariant), and u(t) is a zero- 
mean white noise vector, such that 

Eu(t)u'(r) — KS(t - r) (3) 

where K is an m X m symmetrical matrix and prime de- 
notes the transpose. 

It is desired to calculate the output covariance (an 
autocorrelation, since y(t) is of zero mean) Ey(t)y(r) f 
given by 

(4) 

Let the random variable x(f 0 ), where t 0 is the initial time, 
be statistically independent of u(t). It is well known (Ref. 
5) that the solution of x(f) is given by 


Theorem 1 

A necessary and a sufficient condition that a stationary 
process y(t) be representable as the output of an asymp- 
totically stable, time-invariant, finite, dimensional linear 
model is that its spectral density R(s) be a rational func- 
tion of the form H(s)H(-s), with 


X(t) = #(*, fo)s(fo) + [ $(*, T)B(r)u(T)dr 
Jt a 

where $(£, T ) is the state transition matrix; i.e.. 


r) 

dt~ 


= A(t) *(t, r) 




M(s) 

P( s ) 


( 1 ) 


$(f, t) — I 


( 5 ) 


(?) 

( 7 ) 


for some polynomial 


Substituting x(t) from Eq. (5) into (4) and performing 
some mathematical operations, we obtain (Ref. 3) 


p{$) = s n 4- £ 


i=o 


with all roots in the left half part of the .v-plane and 


M(s) = 

IiQ 


Ey(t) y(r) = C(f)*(f, r) P„(t) C'(r)l(t ~ r) 

+ C(t)P x (tMt, t) C'(f)l(r - t) (8) 

PS) = Ex(tM t) (9) 

where 1(f) is a unit step function. 
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From the dynamical model (Eq. 2), P*(£) can be shown 
to be the solution of the differential equation (Ref. 3) 

P. = AP X + P x A f + BKB f (10) 

where the covariance P*(f 0 ) must be given. 

C. Independence of Estimation Problem of a 
Particular Coordinate System 

In spectral realization, i/(f), given by Eq. (2), is the 
signal without any noise contamination. Often, we receive 
a contaminated observation z(t ), given by 

z(t) = y(t ) + n(t) (11) 

where n(t) is additive noise which is assumed to be un- 
correlated with y(t). Kailath (Ref. 1) and Anderson (Ref. 
6) show that the only information necessary for recursive 
estimation is the knowledge of Ey(t)y(t 4 r) and Ez(t) 
z(t 4- t). That is, the solution of recursive estimation in 
the mean-square sense is independent of the particular 
coordinate system for model z( •) and «/(•) processes; hence, 
a unique solution associated with minimum mean-square 
estimation can be obtained where the models for the 
processes are not given in advance. All these models are 
related to one another by a linear transformation. For 
example, if 

x = Ax(t) 4- Bu(t) 

y = Cx(t) 4 v(t) (12) 

and 

x* = AV(t) + BV(f) 
y = C*x*(t) 4- v(t) (13) 

correspond to the same realization, then there exists a 


linear transformation T(t) such that 


At) = T(t)x(t ) 

(14) 

and 


** = T(tfat) 

(15) 


where x and x* are the estimates corresponding to Eqs. 
(12) and (13), respectively. 


The covariance estimates can be obtained accordingly. 


II. Recursive Image Estimation 

A. Procedure Outline 

The enhancement of images that are characterized only 
by statistical data where the picture contains additive 
noise is considered in this section. The random process 
representing, the output scanner is characterized by the 
output of a dynamical model with white noise input. The 
dynamical model describes the first-order vector Markov 
process. The procedure of Kalman filtering is then utilized 
to recursively determine the minimum mean-square error 
estimate of the image. The result is also extended to ob- 
tain the smoothing of data. Two examples (Ref. 7), one 
with very high SNR, are used to illustrate the effective- 
ness of the procedure. In what follows, the image is 
assumed to be a two-dimensional, stationary correlation 
function of zero mean. Thus, the autocorrelation function 
and the covariance become identical. The statistical in- 
formation about the image and the noise is assumed to 
be known and uncorrelated, and the noise is additive. 

B. Derivation of Autocorrelation Function 
of Scanner Output 

Let us scan a picture horizontally using an optical 
scanner denoted by s(t). Let the horizontal position (a 
continuous variable) be denoted by z, where 0 ^ z ^ Z, 
and the vertical variable by an integer n = 1, 2, •••, N 
representing the nth scanned line. The brightness func- 
tion is defined by b(z y n). Let us assume, without any loss 
of generality that b(z 7 n) is of zero mean. The random 
process b(z , n) is assumed to be wide-sense stationary, 
with the autocorrelation function defined by 

Eb(z 2> n 2 )b(z l7 n x ) = R(z 2 - z l7 n 2 - n x ) = R(z, n) 

(16) 

Assume that the scanner output s(t) has a horizontal 
speed o — l and, without any loss of generality, that the 
vertical movement takes zero time. Let us determine 
Es(t)$(t 4- t) in terms of R(z, n) and Z. The variables t and 
r can be equivalently expressed by 

f “ /T 4- cr, / — 0, 1, 2, — 1, O^cr^T 

r = IT + y, i 1, 0, 1, * • • 

O^t 4 t^NT, O^y^T (17) 
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where T = Z is the time required to traverse one hori- 
zontal line. The scanner output can now be written as 


s(t) = b (<t, / 4- 1), 


/ b (a J r y, i + / *f 1) 

J if <r + y — T 

j fc(«r + y-r,i+/ + 2) 
\ if (j 4- y > T 

(18) 


Now, utilizing Eqs. (16) and (17), we can obtain 


(*) $ (t + r) 


R(y,i) 

if o- + y — T 

R (tx 4- 7 — T y i + / + 2) 
if rx 4~ y > T 


(19) 


It is clear that Es (*) s (t 4- t) is a function of both tr and 
y, or equivalently, of t and t; thus, it must be nonstation- 
ary, The nonstationarity is due to the edge condition. A 
simple check shows that Es ( t ) s(t 4- t) is also periodic 
and a nonseparable function. It can be demonstrated that 
no finite-dimensional linear realization of this nonsepara- 
ble autocorrelation exists (Ref. 3). 


We shall now seek to generate a random process de- 
noted as q (t) such that it has a stationary autocorrelation 
function which approximates the autocorrelation of the 
process s (t). To generate q (t), we proceed as follows. For 
a given t > q ( t ) is defined by 

q(t)=*(jT + f) (20) 

where | is assumed to be uniformly distributed over 
[0, T]. Now we shall prove the following theorem (Ref. 8). 

Theorem 2 

The random process q (£) defined by Eq. (20) is sta- 
tionary. 

Proof 

It is easy to verify that 

Eq (t) = 0 

by the construction of q (t)> 


Next, we must prove that Eq (t) q (t + t) is a function 
of r (or equivalently, y) only. To accomplish this end, we 
calculate the correlation function of the process q [t): 

Eq (t) q(t + t) = E^Es [s (/T 4- |) s (;T 4- £ 4- iT 4- y)] 

= U(fT + i) 

X * (;T 4- <7 + iT 4- y)] dUr (21) 

This equation is obtained by utilizing Eq. (20) and 
t — iT 4- y, which is given by (17), and the fact that | is 
uniformly distributed over [0,T]. The subscripts 5 and £ 
in (21) denote the expectation with respect to s and £, 
respectively. From Eqs. (19) and (21), one obtains 

Eq (t) q (t 4 - t) = R (y, <) d£ 



= ^^R(y,i) +lR(T-y,i + 1) =f(r) 

( 22 ) 

where Eq (t) q (t 4- r) is defined as r (t), which is a func- 
tion of t (or y) only. 

It is interesting to note that the correlation function of 
q (t) } namely r( T ), can also be obtained by averaging the 
autocorrelation function of s (t) over one period. How- 
ever, it is important to mention that such averaging over 
the subintervals of a period may not give rise to a sta- 
tionary autocorrelation function, and furthermore, it may 
not yield an autocorrelation function at all. 

As an example, consider a scalar random process char- 
acterized by a scalar differential equation 

x — — x 4- u 
y (t) = cos (t) x (t) 

where the initial state x (0) — 1/2 and 
Eu ( t ) = 0 

Eu (h) u (f 2 ) = S (t 2 — fi) 

Then, the autocorrelation of x(t) can be obtained as 
follows : 


Ex(f)*(f + r)=-s-exp(— | t|) 
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Thus, Ey (t) y (f 4- r) is given by 

Ey (t) y {t + r) =icos (f)cos(f 4- r)cx p( — |r|) 

which is clearly nonstationary, since the correlation func- 
tion of y (t) depends on both t and t + r and is periodic 
(of periodicity 2ir). However, if we averaged this auto- 
correlation over [0, -tt/4], the resulting average would 
depend on both t and t 4~ r. 

The randomization of £ over the period T has the intui- 
tive appeal that all points of the picture are weighted 
equally. 


resulting in a zero mean sample function. As a first approxi 
mation, let us choose 

R(x,i) = «exp(-^|z| — A%|t|) 

where a, fi h , and are to be determined. Computation 
of the sample power results in a = ft (0, 0) ^ 6.1. The cor- 
relation between two adjacent grid points is calculated as 
5.33, which is a value for R (1/32, 0) or R (0, 1). Hence, 

R(xJ) = 6.1 exp (—4.35 1 3 1 — 0.136|t |) 

The correlation function is obtained by substituting the 
above into Eq. (22), and the plot is shown in Fig. 1. 


The following salient properties of r (t) will be used in 
what follows: 

r(iT) = R(0,t) (23) 

Since R (z, n) is an autocorrelation function, 

R(0,n)^R(z } n) (24) 

Thus, from (22) and (23), 


C. Dynamical Modeling of Image Statistics 

In this section, vve wish to derive a differential equa- 
tion model whose solution has an autocorrelation func- 
tion approximating r(r) given by Eq. (22). Since we 
subsequently intend to utilize a Kalman estimator, we 
seek a dynamical model of the form 

x (t) = Ax (t) + Bu ( t ) 

y(t)=Cx(t) (26) 


r (iT + y) 
r(iT) ~ ’ 


for all i, y 


(25) 


The above properties indicate that, in general, the corre- 
lation function t(r) has a periodic nature. 


Example 1 

Consider a square picture subdivided into a 32 X 32 
grid. Let T = 1 second and v = 1. The signal is a 12 X 12 
square starting at the 13th row and 13th column. Let m 
and n represent specific rows and columns, respectively. 
The above signal is represented by the brightness level 
b(m,n) = 6.1 where the signal exists and —1 otherwise, 


where x(t) is an n-dimensional vector, u(t) is a white 
noise vector, and y ( t ) is the scalar signal whose autocor- 
relation function is f(r). 

The procedure followed is to represent an approxima- 
tion to r (t), denoted by r a (r), as a sum of terms such that 
each term can be easily modeled, since, in general, r( T ) 
may not have a rational bilateral transform. The proper- 
ties of r (t) may be utilized to decompose r (r) into the 
product of two functions h(r) and r(r)//i(r): 

( 27 ) 



Fig. 1. Plot of r(r) and r a (r) (dashed curve) as a function of t 
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where h (t) is chosen to satisfy 


h ( iT ) = R (0, i ); for all i (28) 


Since, in many practical cases, the two-dimensional 
correlation function R(z,i) is a monotonically decreas- 
ing function of *, a natural candidate for h (r) is, in those 
instances, a combination of negative exponentials; i.e., 


ft ( T ) = £ h exp ( — A, 1 t | ) (29) 


The function p (r) is then chosen to be a periodic func- 
tion approximating r(r)/ft(r). The approximate correla- 
tion function is 


r*(r)=h(r)p( r) (30) 

Utilizing Eq s. (23) and (28), it can be seen that the func- 
tion r (t)/^(t) is unity at iT and less than unity for all 
other t; furthermore, from (22) and (29), it is an even 
function. Hence, p(r) is chosen to be an even function 
with period T. Thus, a natural candidate for this func- 
tion is 

V (r) = £ a , cos —L T (31) 

i-o * 

Consequently, an element of the function r a ( r ) has the 
form 

2ni 

liCi exp ( — Xi | T | ) cos ~Y~ T (32) 


and there are (/ 4- 1) I such elements. 


A differential equation model with white noise input 
can be simply constructed (Ref. 8) to model each of these 
terms. Each will be a second-order system except for 
those corresponding to / = 0; i.e., 


which will be of first order. If the white noise forcing 
functions (one being necessary for each i y j pair) are 
chosen to be mutually independent, the collection of all 
these differential equations defines the parameters A, B, C 
and represents the desired model for f a (r). 

In the course of selecting the approximate function 
r a (r), we must choose the coefficients properly, such that 
r a (r) is a correlation function. We shall either guarantee 
that r a ( •) is a positive definite function or equivalently, 
that the spectral density of r a (t) is positive (Ref. 9). 

Example 2 

Using example 1, let us derive a dynamic model for 
r (r). Assume that the desired model has the form given 
by Eq. (26), and further, that 

Eu(t)u(t + rY = K»(r) (33) 

where S (r) is the dirac delta function, the prime denotes 
the transpose, K is a positive definite matrix, and 

(t) y{t + r)^r a ( T ) (34) 

Because of the exponential nature of R(z,i), we choose 

h(r)=R( 0,0)exp(-0.136lr|) (35) 

and 

j 

p(r) ~ 2 cos 2*7 T (36) 

)=Q 

In this example, we use the notation p v instead of 0.136. 

The modeling procedure can be broken down as fol- 
lows. The first term r a (r), namely, 

a 0 exp^-^jrl) 

has tlie bilateral transform 


Uao exp(~-A* |r|) 


(S + fi v ) ( S — /%) 


= HiW 


(37) 
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where the superscript denotes the model corresponding 
to the appropriate term. The coefficients A (2) , B (2) , and 
C <2) are given as 


B (2 > = 


— (2n) 2 + pi -2fi v m 

V 

V2^7„[V (2ir) 2 + rt -2^1 . 

C< 2 > = [1 0] 


In general, the (K + 1) term of ^(t) is a k exp(— /*„|r|) 
cos 2irkT has the bilateral transform given by 


„ , x 2a fc *,[-s 2 + (2fcr) 2 + f4] 

Rk+ ' [S) [(« + *,) 2 + (2fcr) 2 ] [( - * + f**)* + (2&r) 2 ] 


As before, the function Rui(s) can be factored into two 
functions, H*n(s) and Hu i( — s ) : 

V 2g fc ^y[s + V (2krrY + ft 2 ] 

«k + ,(s) - + ^Y + (2 *kY 

V 2ay^[ — s + V (2for) 2 + pi] 


( — s + /JvY + (2fe r) 2 


where 


„ (V V 2a kJ u B [s + V ( 2 M 2 + 

Hut{s) - + ^) 2 + (2kwY 

and the corresponding dynamical model is 

x<*+n = A ( * +1) * ( * +1) (t) + B (fc+1) u ( * +1) (t) 

= C (fc+1> * (ft+1) (t) (40) 

where 


^(»+d = 


— (2fcrr) 2 + Yi - 2 Mv. 


V 2c»m» , 

B <s+1) = i , (42) 

LV2an4.iV (2M 2 + m5-2^]J 


C<*«> = [1 o] 
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It can be seen that the first term of r a (r) is modeled by 
Eq. (38), which is a first-order system, and the subse- 
quent terms by (39), which is the second-order system. 
Thus, to model the (/ + 1) terms of r 0 (r), we need a 


(2/ + 1) -order system. For example, suppose the func- 
tion r a (t) has (/ + 1) terms, then we can incorporate the 
first- and second-order systems into a new system, whose 
parameters A, B, and C are obtained as follows: 


pv 

0 

0 

0 


0 0 

0 1 

~[(2n)* + tf] -2p v 

0 0 


0 


0 


(44) 


0 0 0 
0 0 0 


0 1 

— [( 2 n- (/ + l )) 2 + pi] —2p v 


V 2 o 
0 V 

0 V&h* [V « + A - W] 


o 


0 


0 

0 

0 

0 

V 2 Ojpv 

y[2ajfh, [V (2w/) 2 + - 2fi) _ 


(45) 


C = [1 1 0 1 0 1 0] 


(46) 


Example 3 The second term in the correlation is modeled by 

If in example 2 only three terms of r a (t) are retained, 
i.e., J ~2, the resultant t a (r) can be written as % 2 = x 3 + 0.82 w 2 


Ta (t) = 6.1 exp — 0.136 (r) 22 a i COS 2 irr 

i-o 

If we use the Fourier series for p (t), then a,,, a u and a t 
will be given as 

flo = 0.333; ci\ = 0.405; a 2 = 0.101 


x 3 = — 39.4 x 2 — 0.27 x 3 + 4.92 u 2 

and the third term is modeled in a similar manner. The 
terms u„ u 2> and u 3 represent independent white-noise 
terms, each with zero mean and correlation function 8 (r), 
where 8 is the dirac delta function. The final results are: 


A plot of r a (t) is shown in Fig. 1. The correlation term 
6.1 a 0 exp(— 0.136 |t|) 


is modeled by 


ii = — 0.136 x 2 (t) + 0.732 u. 


A = 


0.136 

0 

0 

0 

0 


0 0 0 

0 1 0 

-39.4 -0.27 0 

0 0 0 

0 0 -157.7 


0 “ 
0 
0 
1 

-0.27 
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" 0.743 0 

0 0.820 

B = 0 4.92 

0 0 

0 0 


C = [1 1 0 1 0] 


0 “ 

0 

0 

0.410 

5.04 


Often, two-dimensional stationary correlation functions 
can be approximated by a combination of two-dimensional 
stationary correlation functions of the form 


To minimize 31 (r), we must minimize 

[r( r)-f.(r)]*dr 

Thus, the minimization of 31 (r) becomes a simple prob- 
lem, and the risk function can be obtained from Ref. 11. 
The procedure is to set the derivatives of 31 (r) with re- 
spect to (ij equal to zero, and the result can be obtained 
as follows (Ref. 11): 


a = a^d 


(52) 


H(x,i) = R( 0, 0) exp ( — At* | x | — | i | 


(47) 


where a is a matrix, whose elements are given by 


Because of the importance of R (; x 9 i) as given by Eq, (47), 
we shall discuss this special autocorrelation function 
below, 



— 2 fjyy | r [ ) COS 2 irkr COS 2 7 rlr dr 


(53) 


Calculating r (r) (given by Eq. 22), one obtains 
r(r) = T T 7 exp(-w|y| -/*,|t|) 


and d is a column vector, whose elements are given by 
*-r T (r) exp ( — I T | ) COS 2-rrfcr dr (54) 


+ y-exp( — | T — y| — p v \i + 1|) (48) Furthermore, the following properties can easily be 

established: 


where 


t = iT + y, 0 — y — T 

Now, let us define a risk function <&(•) such that 


-i: 


£(r)= [r{r)-r a (r)V-dr 


(49) 


J [f(r) ~ f«(r)] z dr = J r 2 (r) - j d (r) dr 

I r 2 (r) dr = lim f r\ (r) dr 
Jo J -4 oo J 0 


(55) 

(56) 


and 

y 2W III. Design of the Estimator 

= z^ ex p(-^'] T i) cos -y‘ r ( 5 °) 

>= ® A. Design of a One-Step Predictor 


We can select the coefficients a, such that the risk func- 
tion 31 (r) is minimized. For simplicity, we shall assume 
that T = 1. It can be shown that 31 (r) can be expressed 
by (Ref. 11) 


Since we intend to utilize a digital computer for the 
estimation process, the model given by Eq. (26) is dis- 
cretized, yielding 


M = 


1 - exp(— 2 pvN) f 1 r , x 

l-«p(-s *.)]. [,(t) 


T a (t)] S dr 

(51) 


*(fc + l) = Ax{k) + Bu(k) 

y (k) = Cx (k) + v ( k ) (57) 
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In addition, the model given by Eq. (57) contains the 
observation noise element v (k), which is assumed to be 
white, with mean zero and variance a 2 . The parameters 
A, B , and C are related to A, B, and C by (Ref. 10) 

A = exp(A^-) 

BKB' = f exp ^A-^-) exp( — As) BKB' 

X exp ( — A's) exp ( A' ~) ds 

C=C (58) 

where K and K are covariances of if (t) and u (7c), respec- 
tively. The sampling interval utilized in the above dis- 
cretization is chosen to be T/N. Thus, there will be N 
observations for each horizontal scan. Since there are 
N horizontal scan lines, the final discrete observation is 
on an N X N grid. 


2(k + 1) = [A - F (k) C] x (k) + F(k)y ( k ) 

P(k + l) = [A-F(k)C]P(k) [A — F(fc)C]' 

+ B KB' + F (k) F' (k) <r 2 
F(k) - AP (k) C' [CP {k) C' + a 2 ]- 1 (59) 

The (one-step predicted) estimate of the image is there- 
fore 

C2(k)±$(k) 

that is, y (fc) is the best estimate of y (k), obtained recur- 
sively in real time, where 1 / ( k ) is the observation associ- 
ated with the grid point immediately ahead of the scanner 
position. 

Example 5 

The signal y ( k ) is generated by using the image de- 
scribed in the preceding example and adding white noise 
with variance <r 2 . Let us define a measure of signal-to- 
noise ratio by 


peak-to-peak variation of signal 

Example 4 ^ & 


Continuing example 3, we obtain 


A — 


BKB' = 


0.996 0 0 0 0 

0 0.983 0.031 0 0 

0 -1.22 0.97 0 0 

0 0 0 0.926 0.03 

0 0 0 -4,77 0.913 

0.02 0 0 0 

0 0.02 0.12 0 

0 0.12 0.60 0 


0 0 
0 0 


0.01 0.07 
0.07 0.49 


C = C = [1 1 0 1 0] 


Utilizing the model given by Eq. (57) with parameters 
given by Eq, (58), a (one-step predictor) recursive esti- 
mator may be designed (Ref, 10). The equations are given 
for the sake of completeness. 


The peak-to-peak variation of the image is 7.1. Two 
values of p are considered here, namely 7,1/3 and 7.1/10; 
the corresponding values of tj ( k ) and their one-step pre- 
dicted values y (k) are shown in Figs. 2a and b and 3a 
and b, respectively. 


B. Implementation of Required Interpolation 

It is clear that image enhancement, from the point of 
view of scanner output, represents an interpolation prob- 
lem; i.e,, it is desired to determine the best estimate of 
y (k), O^k — ZV, given the observation y (0), y (1),- ■ -,t/ (A/). 
In general, the interpolation problem is far more compli- 
cated (Ref. 10) than standard Kalman filtering. However, 
since for the image enhancement considered here the 
length of the data is fixed (N) and, furthermore, the ob- 
servation is usually available for additional repeated 
processing, it is possible to obtain two one-step predicted 
values of y (k), denoted by $ (k) and t/ (k), one by running 
the scanner in one direction starting, for example, at the 
top left corner of the picture and the other by running 
the scanner in the reverse direction starting at bottom 
right corner. Associated with these estimates are estima- 
tion error variances denoted by a 2 (k) — CP (k) C' and 
a 2 (k) — CP (k) C', respectively. The two estimates must 
be combined to yield the optimal interpolated (smoothed) 
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(c) 


(d) 


Fig. 2. Observation and estimates for p = 7/3 


value y*(k). Thus, a brief discussion of combining two 
estimators is warranted. 

Suppose we are given two state estimates, x (t) and 
x(t), of the same state variable x(t). There are two cases 
to consider; either x(t) and lc(t) are correlated or they 
are uncorrelated. We shall combine only the case in 
which both are uncorrelated; i.e., 

E [x — x] [x — x]' — 0 (60) 


In this case, the optimal estimate of x, denoted by x* ( t ), 
is given by 


x* = P*(P‘->£+p-ix) 

(61) 

p* = (P-1 + p-l)-l 

(62) 


where F and F are the error covariances of 'x' and x*, re- 
spectively. Thus, applying _Eq s. (60), (61), and (62) to 
obtain y (k) = Cx and tjf — Cx yields 
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(c) 


(d) 


Fig. 3. Observation and estimates for p = 7/10 


y‘(*) = 


*(*) 

-=nu'J( k ) + 


r(k) + **( k )' 


* ik) •»<*> 


<X“ (Jfc) + 7j 2 (fc) 


(63) 


Example 6 

Considering the preceding example, the covariance 
P (k) in Eq. (59) nearly reaches its steady-state value in 
about two or three scan lines. Consequently, <x(k) ^7r(k) 
for most of the picture, and Eq. (63) reduces to 


y'(k)~j[y(k) + y(k)] (64) 

Equation (64) was implemented, and the results for p = 
7.1/3 and 7.1/10 appear in Figs. 2c and 3c, respectively. 


Careful observation of Figs. 2b and c (or 3b and c) 
reveals a consistent vertical correlation, which is attrib- 
uted to the approximation of r (r) by r a (r) (Fig. 1). The 
effect of this approximation is partially eliminated by 
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transposing the original picture and re-evaluating y'{k)> 
The two values are then averaged and represent Figs, 2d 
and 3d for corresponding values of p. 

IV. Conclusion 

The feasibility of applying recursive (Kalman) filtering 
techniques to image processing has been established. 
Thus, the estimate at any one point does not require 
processing of all the data but only of the information 
stored by the point preceding it. The recursive procedure 


is applicable to those images characterized statistically 
by means and correlation functions. It is important to 
note that the required computational time increases only 
linearly with the size of the picture (number of scan lines, 
number of discrete observations per line). A time-invariant 
dynamical model was chosen, leading to stationary sta- 
tistics for the scanner output. 

Improved performance is expected if the nonstation- 
arity due to the scanner s periodic operation is considered. 
This improvement will be discussed in a later publication. 
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